Old Dominion University Desearch Foundation 


NASA-CR-192980 


DEPARTMENT OF MECHANICAL ENGINEERING AND MECHANICS 
COLLEGE OF ENGINEERING AND TECHNOLOGY 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23529 


GRID SENSITIVITY FOR AERODYNAMIC 
OPTIMIZATION AND FLOW ANALYSIS 


By 




I. Sadrehaghighi, Graduate Research Assistant 
and 

S. N. Tiwari, Principal Investigator 


/bO 3 O 2 — 


Progress Report 

For the period ended December 1992 
Prepared for 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23665 



Under 

Cooperative Agreement NCC1 - 68 
Dr. Robert E. Smith Jr., Technical Monitor 

ACD - Computer Ap plicatio ns Branch 

N93-251] 

Unc 1 as 


(NASA-CR-192980) GRID SENSITIVITY 
FOR AERODYNAMIC OPTIMIZATION AND 
Anril 1993 FLOW ANALYSIS Progress report, 

P period ending Dec. 1992 (Old 

Dominion Univ.) 120 p 


K — =1 

I : 

BS 


G3/02 


0160302 


DEPARTMENT OF MECHANICAL ENGINEERING AND MECHANICS 
COLLEGE OF ENGINEERING AND TECHNOLOGY 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23529 


GRID SENSITIVITY FOR AERODYNAMIC 
OPTIMIZATION AND FLOW ANALYSIS 


By 

I. Sadrehaghighi, Graduate Research Assistant 
and 

S. N. Tiwari, Principal Investigator 


Progress Report 

For the period ended December 1992 
Prepared for 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23665 

Under 

Cooperative Agreement NCC1 - 68 
Dr. Robert E. Smith Jr., Technical Monitor 
ACD - Computer Applications Branch 


April 1993 


FOREWORD 


This is the progress report on the research project “ Numerical Solutions of Three- 
Dimensional Navier-Stokes Equations for Closed-Bluff Bodies” . Within the guidelines 
of the project, special attention was directed toward research activities in the area of 
“Grid Sensitivity for Aerodynamic Optimization and Flow Analysis.” The period of 
performance of this specific research was January 1, 1992 through December 31, 1992. 
This work was supported by the NASA Langley Research Center through Cooperate 
Agreement NCC1-68. The cooperate agreement was monitored by Dr. Robert E. 
Smith Jr. of Analysis and Computation Division (Computer Applications Branch), 
NASA Langley Research Center, Mail Stop 125. 
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ABSTRACT 


GRID SENSITIVITY FOR AERODYNAMIC 
OPTIMIZATION AND FLOW ANALYSIS 

Ideen Sadrehaghighi 
Old Dominion University, 1993 
Director: Dr. Surendra N. Tiwari 

An algorithm is developed to obtain the grid sensitivity with respect to 
design parameters for aerodynamic optimization. Two distinct parameterization pro- 
cedures are developed for investigating the grid sensitivity with respect to design pa- 
rameters of a wing-section as an example. The first procedure is based on traditional 
(physical) relations defining NACA four-digit wing-sections. The second is advocat- 
ing a novel (geometrical) parameterization using spline functions such as NURBS 
(Non-Uniform Rational B-Splines) for defining the wing-section geometry. An inter- 
active algebraic grid generation technique, known as Two-Boundary Grid Generation 
(TBGG) is employed to generate C-type grids around wing-sections. The grid sensi- 
tivity of the domain with respect to design and grid parameters has been obtained by 
direct differentiation of the grid equations. A hybrid approach is proposed for more 
geometrically complex configurations. A comparison of the sensitivity coefficients 
with those obtained using a finite-difference approach is made to verify the feasibil- 
ity of the approach. The aerodynamic sensitivity coefficients are obtained using the 
compressible two-dimensional thin-layer Navier-Stokes equations. An optimization 
package has been introduced into the algorithm in order to optimize the wing-section 
surface using both physical and geometric parameterization. Results demonstrate a 
substantially improved design, particularly in the geometric parameterization case. 
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Chapter 1 


INTRODUCTION 

1.1 Motivation 

Integrated design and optimization of airplane components has become a pri- 
mary objective for most researchers in aerodynamic community. The sudden interest 
can be attributed to the introduction of complex and composite materials required 
by advanced aerospace vehicles, such as National AeroSpace Plane (NASP) and High 
Speed Civil Transport (HSCT) aircraft where, because of extreme flight conditions, 
the interdisciplinary interactions are particularly important. The process requires 
many analyses over a wide range of engineering disciplines. Each analysis is based 
on solving mathematical models describing physical laws associated with a discipline. 
For a vehicle confined to atmospheric flight conditions, the primary engineering dis- 
ciplines are: aerodynamics, structures, control, and propulsion. These disciplines are 
interconnected and affect each other. 

A Multidisciplinary Design Optimization (MDO) would provide the designer 
with sufficient information to predict the influence of a design parameter on all rel- 
evant disciplines. The traditional approach to MDO is to perform the analysis and 
optimization by each discipline in a sequential manner where one discipline uses the in- 
formation from the preceding analysis of another discipline. This tedious and lengthy 
process, although acceptable for loosely-coupled systems, is likely to produce sub- 
optimal results. For systems which are more tightly coupled, a relatively new ap- 
proach is to perform the analysis and optimization at each discipline concurrently. 
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As opposed to a sequential approach, this technique supplies the designer with first 
order (i.e., derivative) information, thus enables him to to predict the influence of 
a design change on all the disciplines involved. The interaction among disciplines 
are achieved by a system of equations known as Global Sensitivity Equations (GSE), 
which communicate the system response due to design perturbations among all dis- 
ciplines. On the local level, within each discipline, the Local Sensitivity Equations 
(LSE) are responsible for similar response. Both GSE and LSE are linear and alge- 
braic, regardless of the mathematical nature of the governing equations associated 
with each discipline. 

A complete design and optimization analysis using all the relevant disci- 
plines is still a formidable task even for an isolated airplane component such as a 
wing or fuselage. The computational cost associated with such analysis can easily 
strain the capabilities of current supercomputers. The magnitude of this problem can 
be best appreciated when a discrete aerodynamic or structural design analysis can 
exhaust the computational capability of a medium size supercomputer. The underly- 
ing problem is the expensive cost of the analysis for each discipline involved. Clearly 
the aerodynamics involve non-linear physics and use of composite materials would 
require non-linear structural analysis as well. For a simple aero-elastic problem, the 
entire system matrix must be simultaneously solved using mostly implicit solvers. 
The extensive computational demand for such coupling of the governing equations, 
will likely limit MDO to only individual components such as a wing or wing-section. 
The cost of optimization operations are relatively small and manageable. Two gen- 
eral directions to overcome these difficulties have been proposed by different research 
groups. The first direction leads toward modifying the existing computational tools in 
order to obtain a relatively cheap and reliable technique for design and optimization. 
The usually favored direct solvers with all their advantages, require extremely large 
computer storage even for 2D applications. 


Recent efforts concentrated on development and implementation of efficient 
iterative techniques and improvement of existing ones. The conditioning of the co- 
efficient matrix, resulting from linearization of the governing equations, are prone 
to affect the convergence criteria and the propagation of error through the system. 
The second direction points to the advent of next generation of supercomputers with 
parallel processing capabilities. Parallel computing would be ideal for MDO analysis 
where each discipline could be assigned to a particular processor for greater efficiency. 
Consequently, the problem formulation and algorithm design (i.e., software develop- 
ment) should change in order to adapt to new computer architecture. Recognizing this 
need, the High Performance Computing and Communication Program (HPCCP) has 
been established by the federal government to confront such problems. This program 
is focused on developing the technology for TeraFLOP computing, an improvement 
of almost 1000 times over current technology. 

For the present, a more realistic task would be to consider a discrete design 
and optimization model for simple configurations. The aerodynamic optimization, 
being an important component of MDO, has become an area of interest for many 
researchers. An essential element in design and optimization of aerodynamic surfaces 
is acquiring the sensitivity of functions of CFD solutions with respect to design pa- 
rameters. Several methods concerning the derivation of sensitivity equations (LSE) 
are currently available. Among the most frequently mentioned are Direct Differ- 
entiation (DD), Adjoint Variable (AV), Symbolic Differentiation (SD), Automatic 
Differentiation (AD), and Finite Difference (FD). Each technique has its own unique 
characteristics. The Direct Differentiation, adopted in this study, has the advantage of 
being exact, due to direct differentiation of governing equations with respect to design 
parameters. The Adjoint Variable, having its roots in structural analysis, produces 
equally impressive results. For Symbolic Differentiation, a symbolic manipulator such 
as MACSYMA can be used to carry out these differentiations. Automatic Differen- 
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tiation, still at preliminary stages, exploits the fact that exact derivatives can be 
computed easily for all elementary arithmetic operations and intrinsic functions. The 
finite difference approach, simple and until recently the most popular, is based upon 
finite difference approximation of the derivatives. In this approach, a design parame- 
ter is perturbed from the nominal value, a new solution is obtained, and the difference 
between the new and old solution is used to obtain the sensitivity derivatives. This 
brute force technique has the disadvantage of being computationally intensive, espe- 
cially when the number of parameters involved is large. 

Design parameters can be classified according to whether or not they are cou- 
pled. Uncoupled design parameters influence the solution independently and would 
be the major contributors to optimization process. These parameters could be geo- 
metric, flow-dependent, or grid-dependent. The geometric design parameters specify 
the primary shape of a typical aerodynamic surface. Flow-dependent parameters are 
usually free-stream conditions such as free-stream Mach number or angle of attack. 
The grid-dependent parameters, relatively new in aerodynamic optimization, affect 
the interior and boundary grids; therefore, influencing the solution and optimiza- 
tion process. Traditionally, geometric parameters are considered the most affluent 
in aerodynamic optimization, although, optimization with respect to other design 
parameters is gaining respectability. For optimization with respect to geometric de- 
sign parameters, a perturbation in parameters affect the surface grid and the field 
grid which, in turn, affect the flow-field solution. There are two basic components in 
obtaining aerodynamic sensitivity. They are: (1) obtaining the sensitivity of the gov- 
erning equations with respect to the state variables, and (2) obtaining the sensitivity 
of the grid with respect to the design parameters. The sensitivity of the state vari- 
ables with respect to the design parameters are described by a set of linear-algebraic 
relation. These systems of equations can be solved directly by a LU decomposition 
of the coefficient matrix. This direct inversion procedure becomes extremely expen- 
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sive as the problem dimension increases. A hybrid approach of an efficient banded 
matrix solver with influence of off-diagonal elements iterated can be implemented to 
overcome this difficulty. 

1.2 Literature Survey 

The literature on sensitivity analysis and optimization is quite extensive. 
The pioneering work on sensitivity analysis for MDO started with a plea from So- 
bieski [1,2]* to the CFD community for extending their present capabilities to include 
sensitivity analysis of aerodynamic forces. Yates [3] developed an analytical approach 
using an implicit differentiation in combination with linearized lifting-surface theory 
to evaluate the sensitivity coefficients. This can be used as a benchmark criteria 
for assessing the accuracy of approximate methods. A semi-analytical technique, us- 
ing linear unsteady aerodynamics, has been applied to an isolated wing-section and 
rotating propfan blades by Mur thy and Kaza [4]. Some aeroelastic analysis for a 
transport wing has been investigated by Grossman et al. [5], where a coupled aero- 
dynamic and structure model influences the design. Livne et al, [6] and a few other 
researchers focus on more complex interactions such as inclusion of active controls on 
the overall optimization process. Elbanna and Carlson [7] developed a quasi-analytical 
technique for evaluating wing-section aerodynamic sensitivity coefficients in transonic 
and supersonic flight regimes. Later, they extended the technique to 3D full potential 
equations using the symbolic manipulator MACSYMA to obtain the sensitivity coeffi- 
cients. The procedure was applied to a ONERA M6 wing planform with NACA 1406 
wing-sections [8]. For non-linear aerodynamics, most of the efforts are concentrated 
on involvement of CFD for both flow and sensitivity analyses. Baysal and Eleshaky 
[9,10] presented an aerodynamic design strategy using direct differentiation of Euler 
equations. The procedure was applied to design a scramjet- afterbody configuration 
lumbers in brackets indicate references 
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for an optimized axial thrust. This scheme was later extended to include domain de- 
composition capabilities in order to reduce the computational costs associated with 
complex configurations [11]. 

Taylor et al. [12] conducted a feasibility study of sensitivity analysis in- 
volving Euler equations. The method was successfully applied to two test problems, 
including a subsonic nozzle and a supersonic inlet. The sensitivity derivatives are 
obtained by direct differentiation of governing equations with respect to geometric 
design parameters. The authors later expanded the formulation to include thin-layer 
Navier-Stokese equations and a optimizer. Aerodynamic sensitivity derivatives were 
obtained for an internal flow through a double-throat nozzle and an external flow 
over a NACA four-digit wing-section [13]. Both geometries were optimized with new 
design having a significantly improved performance. The flow and flow sensitivity 
analysis module (ANSERS), developed by these authors, have been implemented in 
this study. Burgreen et al. [14] improved the efficiency of an aerodynamic shape op- 
timization on two fronts. The first improvement involves replacing a previously grid 
point-based approach for surface parameterization with a Bezier-Bernstein polyno- 
mial parameterization. The second improvement includes the use of Newton’s method 
instead of familiar and expensive Alternating Direction Implicit (ADI) technique to 
calculate the flow solutions. Other notable schemes include variable- complexity de- 
sign strategies, developed by Hutchison et al. [15,16], to combine conceptional and 
preliminary-design approaches. The strategy has been used to optimize the HSCT 
wing configuration. Verhoff et al. [17,18] developed a method for optimal aerody- 
namic design of wing-sections using analytically computed aerodynamic sensitivities. 
The scheme also utilizes Chebyshev polynomials together with parametric stretching 
functions to define camber and thickness distribution of wing-section. Due to analyt- 
ical parameterization of surface, the package produces an efficient optimal results. 
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1.3 Objectives of Present Study 

After reviewing relevant literature, it is apparent that one aspect of aerody- 
namic sensitivity analysis, namely grid sensitivity, has not been investigated exten- 
sively. The grid sensitivity algorithms in most of these studies are based on struc- 
tural design models. Such models, although sufficient for preliminary or conceptional 
design, are not acceptable for detailed design analysis. Careless grid sensitivity eval- 
uations, would introduce gradient errors within the sensitivity module, therefore, in- 
fecting the overall optimization process [19]. Development of an efficient and reliable 
grid sensitivity module with special emphasis on aerodynamic applications appears 
essential. 

Unlike aerodynamic considerations, the grid sensitivity analysis has been 
used on structural design models for a number of years. In this context, grid sen- 
sitivity can be thought as perturbation of structural loads, such as displacement or 
natural frequency, with respect to finite element grid point locations [20]. Two basic 
approaches have been cited for grid sensitivity derivatives. The first approach, known 
as implicit differentiation, is based on implicit differentiation of discretized finite ele- 
ment system. The other, which is based on the variation of continuum equations, is 
known as variational or material derivative approach. The main objective here is to 
develop a fast and inexpensive method for grid sensitivity to be used on an automated 
aerodynamic optimization cycle. 

Among two major classes of grid generation systems (Algebraic, Differen- 
tial), algebraic grid generation systems are ideally suited for achieving this objective. 
The explicit formulation, resulting in a fast and suitable grid, enables direct differen- 
tiation of grid coordinates with respect to design parameters [21,22]. The underlying 
effort here is to avoid the time consuming and costly numerical differentiation. In 
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addition, the analytical derivatives are exact, a desirable feature for sensitivity anal- 
ysis. An important ingredient of grid sensitivity is the surface parameterization. The 
most general parameterization would be to specify every grid point on the surface 
as a design parameter. This, although convenient, is unacceptable due to high com- 
putational cost. It is essential to keep the number of parameters as low as possible 
to avoid a surge on computational expenses. An analytical parameterization, may 
alleviate that problem but it suffers from lack of generality. A compromise would be 
using spline functions such as a Bezier or Non-Uniform Rational B-Spline (NURBS) 
function to represent the surface [23,24]. In this manner, most aerodynamically in- 
clined surfaces can be represented with few control (design) parameters. 

Another important aspect of grid sensitivity, grid sensitivity with respect 
to grid parameters, also deserves more attention. This concept, leading to grid opti- 
mization, can be used to enhance the quality of grid in optimization cycle, resulting 
in better flow analysis and faster convergence rates. The accuracy, stability, and gen- 
eral reliability of most flow solvers for any problem may be strongly influenced by 
the choice and quality of grid. This implies that the problem of generating a suit- 
able grid is no longer only restricted to generating a valid grid, but also the related 
issue of manipulating that grid to achieve certain objectives. Among those objec- 
tives, grid smoothness, orthogonality, clustering, and far-field boundary location are 
considered most significant. For example, grid smoothness is an important property 
since an abrupt change in grid size may prompt inaccuracy, ill-conditioning, and in- 
stability in the flow solution. The orthogonality factor can play an essential role in 
finite-difference computations where near orthogonality of grid lines are desirable. 
Also, the accuracy and efficiency of most computational schemes are enhanced by 
grid clustering in regions of high gradients (e.g., boundary layer, shocks, etc.). The 
far-field boundary location has been identified as a dominant factor in influencing the 
solution for a fixed initial conditions. For example, previous investigations indicate 


that for a symmetrical wing-section, the error in lift coefficient has an inverse radial 
dependency on the boundary extent [33]. As required by most optimization tech- 
niques, the sensitivity of the grid with respect to those parameters influencing these 
objectives is required. 

The organization of this study is as follows. The physical and geometric 
representations of a typical model are derived in Chap. 2. The grid generation algo- 
rithm and boundary grid distribution are developed in Chap. 3. Chapter 4 discusses 
the theoretical formulation and aerodynamic sensitivity equation. The method of 
solution is provided in Chap. 5. The results are presented and discussed in Chap. 6. 
Finally, some concluding remarks are provided in Chap. 7. 


Chapter 2 


PHYSICAL MODEL 

2.1 Wing-Section Example 

The physical model considered for this study is an isolated wing-section 
since much research has been devoted to its development and representation. This 
design is essential for the performance of an advanced aircraft for both subsonic and 
supersonic speeds. Other applications could be helicopter rotor blades, and high 
performance fans. Two approaches have been chosen to generate the desired wing- 
sections. The first approach is a physical (i.e., analytical) representation resulting in 
classical NACA four-digit wing- sections. The second approach is a geometric (i.e., 
approximative) representation of the wing-sections using NURBS. 

2.1.1 Physical Representation 

The NACA four-digit wing sections are examined for grid-generation param- 
eterization. Families of wing sections are described by combining a mean line and a 
thickness distribution. The resultant expressions possess the necessary features that 
suit the problem, mainly the concise description of a wing section in terms of several 
design parameters. Reference 25 provides the general equations which define a mean 
line and a thickness distribution about the mean line. The design parameters are: 
T = the maximum thickness, M = the maximum ordinate of the mean line or cam- 
ber, and C = chordwise position of maximum ordinate. The numbering system for 
NACA four-digit wing-section is based on the geometry of the section. The first and 
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second integers represent M and C respectively, while the third and fourth integers 
represent T. Symmetrical sections are designated by zeros for the first two integers, 
as in the case of NACA 0012 wing-section. Figure 2.1 provides a schematic of the sec- 
tion definition. The ^-coordinate is mapped into the chord line x = x(r) = x(/i(0) 
covering both the top and bottom of the section. Details of mapping will be discussed 
in the next chapter. The mean line equation is 


M 

*(») = gj(*7* - * ), x<C 

/ \ _ ~ ~f~ ^^ >a ' ~ a ' 2 ) ■> fl 

y c (x) — M (i _ > — 


The section thickness is given by 


y T (x) = — (0.2969x* - 0.126X - 0.3516x 2 + 0.2843x 3 - 0.1015x 4 ). 
0.2 


( 2 . 1 ) 

( 2 . 2 ) 


(2.3) 


The section coordinates are 


zi(i'.Pi) = 1 !/i(r, P,) = *«(*) ± yr(x) 


(2.4) 


where P\ represents the vector of independent parameters to be defined later. 
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2.1.2 Geometric Representation 

Another approach for representing a wing-section model is using a spline 
function to approximate the surface. The most commonly used approximative repre- 
sentation is the Non-Uniform Rational B-Spline (NURBS) function. The NURBS pro- 
vide a powerful geometric tool for representing both analytic shapes (conics, quadrics, 
surfaces of revolution, etc.) and free-form surfaces [26-28]. The relation for a NURBS 
curve is 


X(r) 


mr=o ^*,p( r ) u; « 




(2.5) 


i = 0, ,n 

where X(r) is the vector valued surface coordinate in the r-direction, D,- are the con- 
trol points (forming a control polygon), u>i are weights, and Ni tP {r) are the p-th degree 
B-Spline basis function defined recursively as 


NM = 


N, 

r — r 


■ , , f 1 r,- < r < r,+, 1 

«,°v; ^ q otherwise ) 


-Wi,p-l(r) + _ r<W AT, tlJ -,(r). 


r i+p - r,- ’ n+p+i - r i+i 

The r< are the so-called knots forming a uniform knot vector 


*} 


( 2 . 6 ) 


(2.7) 


r — < o • • • q , r p ^-i , ...., r m -p_i, /> • • • b ^ 
l p+i p+i 

where the end knots a and b are repeated with multiplicity p + 1. The degree, p, 
number of knots, m T 1 , and number of control points, n + 1, are related by 


m = n + p + 1. 


(2.8) 


i, ill mill lilli iui*iiJi| 8IM IM.iiil '! ' ■ »< il t • m 1 m,xiM i nkih i< 
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For most practical applications the knot vector is normalized and the basis function 
is defined on the interval (a = 0, b = 1). Equation (2.5) can be rewritten as 


x(r) ~S ft "’ (r)Di Ri ' p{r) nJ'wjwu.,- 


i = 0, 


(2.9) 


where Ri, p (r) are the Rational Basis Functions , satisfying the the following properties 
among many others found in [22] 


= i (2.io) 

t=0 

Three options are available to define a wing-section using the NURBS al- 
gorithm. In the first option, the camber line is defined by a NURBS curve using 
three control points. The thickness distribution, Eq.(2.3), is then added and sub- 
tracted to the camber. The first and last control points are fixed for the section 
chord. The design parameters using this option are the location of the middle control 
point, its weight, and the maximum thickness T as shown in Fig. 2.2. Figure 2.3 
shows the corresponding quadratic basis function (p=2,n=2) with weights set to 1 
(i.e., u = l,i = 0,2). The choice of number of control points is a trade-off between 
complexity and functionality [4]. Figures 2.4 and 2.5 illustrate the effect of increasing 
the number of control points on camber, wing-section, and the basis function. 

The second option is to define both camber and thickness distribution curves 
using NURBS representation. The new wing-section can be obtained using Eq. (2.4). 
Both camber and thickness distribution curves are defined using three control points. 
This approach, although promises more design control, it also increases the number 
of design parameters as shown in Fig. 2.6. 

The third option is to bypass the camber and thickness distributions com- 
pletely and control the wing-section directly with NURBS control points and weights. 
Figure 2.7 illustrates a seven control point representation of a wing-section. The 
points at the leading and trailing edges are fixed. Two control points at the 0% chord 
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are used to affect the bluntness of the section. The movement of control points as 
shown in Fig 2.8 creates the effect of camber in the wing-section. The cubic basis 
function (p=3) using this approach is shown in Fig. 2.9 with weights set to 1. 
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(c) Schematic of wing-section 

Fig. 2.2 Wing-section specification using NURBS (option 1). 
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Fig. 2.4 Effects of increasing the number of control points on camber and wing-section 
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Fig. 2.5 Effects of increasing the number of control points on a quadratic basis function 
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(a) Thickness distribution 


■ Control points 

— Control polygon 

— Camber line 






Fig. 2.7 Seven control points wing-section specification using NURBS (option 3) 



Fig. 2.8 Effect of control point movement. 


21 





pjjlfF* n 


Chapter 3 


GRID GENERATION 

3.1 Introduction 

In order to study the flow-field around any aerodynamic configuration, a sys- 
tem of nonlinear partial differential equations must be solved over a highly complex 
geometry [29]. The domain of interest should be descretized into a set of points where 
an implied rule specifies the connectivity of the points. This discretization, known 
as grid generation, is constrained by underlying physics, surface geometry, and the 
topology of the region where the solution is desired [30-32]. A poorly constructed 
grid with respect to any of the above constraints, may fail to reveal critical aspects 
of the true solution. 

The discretization of the field requires some organization in order for the 
solution to be efficient. The logistic structure of the data such as grid spacing, the 
location of outer boundaries, and the orthogonality can influence the nature of the 
solution [33,34]. Furthermore, the discretization must conform to the boundaries of 
the region in such a way that boundary condition can be accurately represented .[35]. 
This organization can be provided by a curvilinear coordinate system where the need 
for alignment with the boundary is reflected in routine choice of cartesian coordinate 
system for rectangular region, cylindrical coordinate for circular region, etc. This 
curvilinear coordinate system covers the field and has coordinate lines coincident 
with all boundaries. To minimize the number of grid points required for a desired 
accuracy, the grid spacing should be smooth, with concentration in regions of high 
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solution gradients. These regions may be the result of geometry (large surface slopes 
or corners), compressibility (entropy and shock layers), and viscosity (boundary and 
shear layers). A complex flow may contain a variety of such regions of various length 
scales, and often of unknown location. 

Two primary categories for arbitrary coordinate generation have been iden- 
tified. There are algebraic systems and partial differential systems. The algebraic sys- 
tems are mainly composed of interpolative schemes such as TVansfinite Interpolation 
[36], Multi-Surface Interpolation [37], and Two-Boundary Interpolation techniques 
[38]. The basic mathematical structure of these methods are based on interpola- 
tion of the field values from the boundary. For partial differential equation systems, 
a set of partial differential equations must be solved to obtain the field values. The 
differential methods may be elliptic , parabolic, or hyperbolic, depending on the bound- 
ary specification of the problem. Each of these grid generation systems has its own 
advantages and drawbacks depending on geometry and application of the problem. 
Algebraic generating systems offer speed and simplicity while providing an explicit 
control of the physical grid shape and grid spacing. However, they might produce 
skewed grids for boundaries with strong curvature or slope discontinuity. Partial dif- 
ferential systems, although offer relatively smooth grids for most applications, are 
computer intensive, specially for three-dimensional cases. An alternative, a common 
practice in recent years, has been to originate the grid using an algebraic system and 
then smooth the field using a differential system. Such hybrid approach proven to be 
successful and cost effective for most applications. 

An array of general purpose grid generation softwares have been emerged 
over past few years. Among many others, the GRAPE2D of Sorenson [39], the EA- 
GLE of Thompson [40], and GRIDGEN by Steinbrenner et al. [41], are the most 
widely used. The GRAPE2D solves Poisson’s equation in two-dimension and uti- 
lizes a novel approach for determination of the boundary control functions. The 
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EAGLE code combines techniques in surface grid generation as well as two or three- 
dimensional field grid generation. The GRIDGEN series is a more recent appearance 
with both algebraic and differential generation capabilities on an interactive environ- 
ment. Another new arrival, called ICEM/CFD, has the capability of combining a 
full Computer Aided Design system (CAD), with grid generation module [42]. This, 
provides an efficient and also quick procedure to reflect the CAD model changes on 
grids. Most of these packages furnish a host of options with a high degree of flexibility. 
However, intelligent use of the majority of these options requires the user to be well 
versed in current grid generation techniques. 

Due to directness and relative simplicity of algebraic systems, the remain- 
der of this chapter would be devoted to their development. The relevant aspects of 
algebraic generation system such as boundary coordinate transformation, mapping, 
boundary discretization, and surface grid generation are discussed in the following 
sections. 


3.2 Boundary-Fitted Coordinate Transformation 

Structured algebraic grid generation techniques can be thought of as trans- 
formation from a rectangular computational domain to an arbitrarily -shaped physical 
domain as shown in Fig. 3.1 [43]. The transformation is governed by vector of control 
parameters, P, and can be expressed as 


where 


x(e,*.c, p) 


r*(*,7.c.p)i 

y(£,»bC> p ) ► 

U^,C,P)J 


Q<(<1, 0 < »7 < 1, and 0 < < < 1. 


(3.1) 


The control parameter P, is composed of parameters which control the primary shape 
of the boundary (design parameters), and parameters which control the grid 
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(grid parameters). A discrete subset of the vector-valued function X-(£i,*lj,(k, P) = 
X{x y z }Jj k = X* is a structured grid for & = = 77 ^ » where 

t = 1,2,3* j = 1,2,3, ••• , M and A: = 1,2,3, • • • , N. 

The orientation of the computational coordinates relative to physical coordi- 
nate, known as grid topology , is an important aspect of the transformation procedure. 
In order to establish a grid topology for any geometry, it is essential to examine 
each component separately [44,45]. For any given geometry, there are several pos- 
sible topologies with different characteristics in terms of efficiency, coordinate cuts, 
singularities, etc. For example a typical wing-section geometry, may have at least 
three types of different topologies (e.g., C-, 0-, or H-types) as illustrated in Fig. 3.2. 
The C- and O-type topologies usually produce the most efficient grid. For present 
study a C-type topology has been chosen and the mapping is shown in Fig. 3.2(a). 
This topology produces no singularity and it is relatively simple to implement. For 
wing-sections with sharp noses, a H-type topology would be more appropriate. For 
more complex geometries, selection of different computational coordinate systems for 
different regions of physical domain might be required. In this case, physical domain 
is mapped into several computational sub-domains, where each sub-domain is refe'red 
as a block. Therefore, it is possible to have a boundary-fitted coordinate system for 
a highly complex configurations. For example, a typical airplane geometry has two 
main components: the fuselage and the wing. A fuselage has a circular like cross- 
section which suggests a natural O-type (cylindrical coordinate) grids. This topology 
produces a nearly orthogonal grid with one line polar singularity at the nose. For the 
streamwise direction, it is feasible to have either a C-type or a H-type grid depending 
on the slope of the nose. For a fuselage with small nose slope, a H-type grid in the 
streamwise direction would be more appropriate. A wing has its own natural coordi- 
nates which are usually not compatible with the fuselage’s coordinate system. It is 
possible to generate a H-, 0-, or a C-type grids in the streamwise direction, and a C- 
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or a H-type in crosswise direction. To maintain a minimum of C° continuity at the 
interfaces, it is essential to select a compatible topology for the wing and fuselage. 
For most cases it is conceivable to generate a single block grid about these compo- 
nents, but this grid tends to be skewed for any practical purposes. A dual-block grid 
possesses much less skewness than a single-block grid. It consists of two large blocks, 
one covering the top portion of the physical domain, and the other covering the bot- 
tom portion of the physical domain. The dual-block topology is a direct consequence 
of using a H-type grid for the wing of zero wing-tip area. Figure 3.3 illustrates the 
mapping of a generic airplane geometry using a dual-block topology. A C-0 type grid 
have been chosen for a fuselage while the wing, horizontal, and vertical tails mapped 
' to a H-H type grids. 

Once the grid topology has been selected, then the grid on the boundary can 
be generated. The boundary discretization techniques will be discussed in following 
section. 


3.3 Boundary Discretization 

Before generating the interior grid, the grid-point distribution along the 
boundary edges should be computed. A discrete uniform distribution of the com- 
putational coordinate can be mapped into an arbitrary distribution of the physical 
coordinate, using a specified control function. The essence of mapping is that the 
abscissa corresponds to the percentage of grid points and the ordinate corresponds 
to a particular control function which, in turn, relates to the geometric definition of 
the physical domain. The control functions must be monotonic in parametric space, 
and can be computed by an analytical function or by a numerical approximation. 
Analytical functions are generally limited to simple boundary curves. However, a 
complex boundary can be decomposed to several sections, and analytical functions 


Fuselage 
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can be used for each section [46]. The advantage of analytical functions are their 
simplicity and the guaranteed monotonicity. Exponentials are the most widely used 
analytical functions and can be expressed as 

B,_,( v'f'jj-J, I _ j 

r(0 - Y,. i + [Yj - Y,. ,1 eB ,_ t _ j for X,-i<(<X, (3.2) 

where 

0 < X i} Yi < 1 , 0 < (,r(0 < 1 and t = 2,....,m. 


The integer m represents number of control points with coordinates The 

quantity £,_i , called the stretching parameter, is responsible for grid density. Speci- 
fying B\, values of Z?,> 2 are obtained by matching the slopes at control points. This, 
guaranteeing a smooth grid transition between each region, can be accomplished us- 
ing Newton’s iterative scheme which is quadratically convergent. 

The exponential function, while reasonable, is not the best choice when the 
variation in grid spacing is large [35,47]. The truncation error associated with the 
metric coefficients is proportional to the rate of change of grid spacing. A large vari- 
ation in grid spacing, such as the one resulting from exponential function, would 
increase the the truncation error, hence, attributing to the solution inaccuracies [35]. 
A suggested alternative to exponential function has been the usage of hyperbolic sine 
function given as 


rtf) = K-i + [Y - x] 


sinh(Bi-i) 


for Xi-x <i<Xi 


(3.3) 


where 

0 < Xi,Y < 1 , 0<{,r(O<l and t = 2,....,m. 


The hyperbolic sine function gives a more uniform distribution in the immediate 
vicinity of the boundary, resulting in less truncation error. This criteria makes the 
hyperbolic sine function an excellent candidate for boundary layer type flows. A more 
appropriate function for flows with both viscous and non-viscous effect would be the 
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usage of a hyperbolic tangent function such as 

for Xi-i < { < Xi 
(3.4) 

where 

0<Xi,Yi<l , 0<£,r(£)<l and * = 2,....,m. 

The hyperbolic tangent gives more uniform distribution on the inside as well 
as on the outside of the boundary layer to capture the non-viscous effects of the so- 
lution. Such overall improvement, makes the hyperbolic tangent a prime candidate 
for grid point distribution in viscous flow analysis. Figure 3.4 illustrates these distri- 
bution functions and the corresponding grids on a unit square. 

Similarly, a numerical approximation can be used to compute the grid-point 
distribution on a boundary curve. This approach is widely used for complex con- 
figurations and care must be taken to insure monotonicity of the distribution. For 
example, the natural cubic spline is C 2 continuous and offers great flexibility in grid 
spacing control. However, some oscillations can be inadvertently introduced into the 
control function. The problem can be avoided by using a smoothing cubic spline tech- 
nique and specifying the amount of smoothing as well as control points [46]. Another 
choice would be the usage of a lower order polynomial such as Monotonic Rational 
Quadratic Spline (MRQS) which is always monotonic and smooth [29]. Other advan- 
tage of MRQS over cubic spline is that it is an explicit scheme and does not require 
any matrix inversion. 


r(0 = 1 + [Yi- Yi . ,] { 1 + 


i- 1] j] 


- !)) 


tanh^f 1 ) 




3.4 Transfinite Interpolation 

The dominant algebraic approach for grid generation is the Transfinite In- 
terpolation scheme. The general methodology was first described by Gordon [36] , 
and 
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there have been numerous variations applied to particular problems. The method- 
ology can be presented as recursive formulas composed of univariate interpolations 
[45] or as the Boolean sum of univariate interpolations [43]. Here, we follow the 
Boolean sum representation; but, for brevity, we restrict the process to two dimen- 
sions and omit some of the details that can be found in [43]. Also, to be consistent 
with the example considered below, the parameterization is restricted to functions 
and first derivatives at the boundaries (£ = 0,1) and (17 = 0,1) and control in the 
interpolation functions. The transformation is 


X(£, ?7,P) = U©V = U + V — UV (3.5) 


where 


u = EE“?(f.Po) 


7=1 n= 0 


yxt&.T.p;) 

d(" 


(3.6) 


and 

2 1 Qm 

v=EEW- 

J= 1 m=0 

The term UV is the tenser product of the two univariate interpolations and can be 
expressed as 


d Tf m 


(3.7) 




UV = EEE£ «?K.pSM?(i.p 5 )- — 

M J=ln=0m=0 ' 


(3.8) 


The boundary curves and their derivatives and - — y J = 1,2 m, 

n = 0, 1) are blended into the interior of the physical domain by the interpolation 
functions a/(£,Po) and /3j(»?,P2). The boundary grid, the derivatives at the bound- 


ary grid and the spacing between points are governed by the parameters {Pj P/j . 

T 

The interpolation functions are controlled with the parameters |Pq Pq} . The en- 
tire set of control parameters can be thought of as a vector |Pq Po Pj P/} = P* 

There are numerous algebraic grid generation techniques which can be de- 


duced from transfinite interpolation formulation. The most successful techniques, 
however, have been those that provide adequate othogonality control and grid spacing 
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control with reasonable function complexity [48]. A prime example of such technique 
would be the two-boundary technique described by Smith [46], using Hermite cubic 
interpolation functions in one coordinate direction between two opposing boundaries. 
In the remaining coordinate direction, linear interpolation between opposing bound- 
aries is specified. This is the technique employed in this study and its development 
and implementation will be the subject of next section. Another notable technique is 
the application of transfinite interpolation using Lagrangian interpolation functions, 
where two, three, or four surfaces are specified in each direction for better grid control 
[29,49]. Detailed formulation of this technique applied to generic airplane configura- 
tion of Fig. 3.3a is outlined in Appendix A. The multi-surface method of Eiseman [48] 
is also another popular grid generation tool. It is a very flexible univariate scheme 
and is similar to Bezier and B-Spline approximation, where the parameters defining 
the surface are not necessarily on the surface itself [50,51]. 

3.5 Two- Boundary Grid Generation Technique 

An interactive univariate version of Eq.(3.5) using only the normal inter- 
polant V is developed. This, called Two-Boundary Grid Generation Technique 
(TBGG), matches both the function and its first derivative at two opposing bound- 
aries. An analytical approximation of the physical coordinates can be expressed as 

x = x,(r,P«)/9f((,PJ) + n(r) d * ,( £ Pl) /)!(t,P!) 

+ x a (x, PM(t, PS) + W , PS) (3.9) 

v = Mr,Pi)0?(i,PS) + 

+ PS)#((, PS) + SM ^W t-PS) (3.10) 


where 


/?i (f j Pq) = 2f 3 — 3f 2 + 1, 
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f3l(t,P v 0 ) = t 3 -2t 7 + t, 

#(*, PS) = -2< 3 + 3< 2 , 

and 

0 < t < 1. 

Figure 3.5 presents the parametric representation of the boundaries and the 
cubic connecting function of Eqs.(3.9) and (3.10). The Hermite cubic blending func- 
tions are shown in Fig. 3.6. 

Five functions r = /i(£)> s = /*(£), R(r ) = Kif 3 (£) f S(s ) = KifiiQ, 
t = f 6 (r/), and their implied defining parameters control the grid spacing on the 
boundaries and the interior grid. Functions r and s define the grid spacing for lower 
and upper boundaries respectively, while R(r) and .S'(s) specify the orthogonality 
along those boundaries. The parameter t defines the grid distribution for the con- 
necting curves between the two boundaries. The quantities K\ and K? are parameters 
that scale the magnitude of the orthogonality at the boundaries. Increasing K\ and 
I <2 extends the orthogonality of the grid into the interior domain. Excessively large 
values of K\ and I(i can also cause the grid lines to intersect themselves, which is 
not a desirable phenomenon. 

For a wing-section example, the domain has been decomposed into an up- 
per and lower section as shown in Fig. 3.7. For the streamwise direction, a bi- 
exponential distribution function with Bi = 4.5 and inflection point coordinates 
{X = 0.3, Y = 0.2} has been chosen for wing-section surface. The relationship, using 
Eq. (3.2) for m = 3, can be expressed as 

r e Bi(J}) Ji 

r tipper = riower = Y | gB, _ j — j 0 < ^ < X 
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r e B *(b=x) — n 

r upper = r i ower = Y + (1 - Y) [ j X < £ < 1. (3.11) 

The stretching parameter, # 2 , is obtained by matching the slopes of the two exponen- 
tial functions at {X, F}. For NACA four-digit wing-section , the section coordinate 
{xi(r, Pi),yi(r, Pi)}, is acquired by invoking Eqs. (2.1-2.4) with 

X — Tapper = Flower' (3.12) 


For a geometric NURBS representation, Eq. (2.9) can be used with 


T — T upper — flower 


(3.13) 


in conjunction with a pre-determined set of control points, D,-, and corresponding 
weights, cu,-. At the wake region, a single-exponential distribution function (m = 2) 
is chosen and, again, the grid continuity is preserved with matching the slopes of the 
distribution functions at the interface. Figure 3.8 shows the resultant discretization 
for upper portion of the domain. For the normal direction, a single hyperbolic tangent 
function such as 

tanh[^(ri — 1 )] 


t = 1 + 


0 < T} < 1 


(3.14) 


tanf i(^) 

with i ?3 = 3.0 is used, concentrating the grid near the wall for capturing the boundary 
layer effects. The grid orthogonality at surface for this particular example is obtained 
using the components of unit normal vector 


dxi (r, P^) 


= ^ sinO 


dj/i(r,Pj) 


= ±CO30 


0 = tan 


■i( dyi(r.Pi) 

\d*,(r,Pj). 


(3.15) 


8t St ' \0*,(r,P5), 

with a constant scaling parameter, (i.e., R(r) = K\ = 2.0). For far-field boundary 

{x 2 (s,P 2 ),j/ 2 ('S,P 2 )}> essentially the same distribution as surface boundary can be 
adopted. The orthogonality at the far-held boundary has been ignored (5(s) = = 

0.0). Figures 3.9 and 3.10 show sample grids for NACA 0012 and NACA 8512 sections 
using this procedure. 
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Fig. 3.9 Example grid for NACA 0012 wing-section 



Fig. 3.10 Example grid for NACA 8512 wing-section 


Chapter 4 


THEORETICAL FORMULATION 


4.1 Generic Sensitivity Equation 

An implicit representation of a physical system can be modeled mathemat- 
ically as 

F(H,G(H)) = 0 (4.1) 

where G and H are dependent and independent variables, respectively. The function 
F can have algebraic, differential, integral or integral- differential characteristics. 
The quantities G and H can be either scalar or vector depending on the nature of the 
physical model. The sensitivity of G with respect to H can be obtained by implicit 
differentiation of Eq. (4.1) 


f 5F 1 _ JdF] I dG\ 
\dH J [dG\\dH) 


(4.2) 


The coefficients, {fjj} and [|^], can be obtained, provided that the solution to 
Eq.(4.1) is known. Equation (4.2), now a set of algebraic equations, can be eas- 
ily solved for the sensitivity derivative, {fjj}- If {Is} an( I [§g] are n °f ava dable, a 
finite difference approach can be adopted. The central difference approximation of 
{ilf} can Revised as 


r dG_ 

\dH 


G(H + AH) - G(H - AH) 

2AH 


(4.3) 


where AH is a small perturbation of a specified parameter. Although the implemen- 
tation of the finite difference approach is comparatively easy, it has the disadvantage 
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of being computationally expensive. Also, the choice of AH is crucial for accuracy 
of the derivative. A large values of AH may lead to inaccurate derivatives while a 
small value may result in round-off errors. 


4.2 Aerodynamic Sensitivity Equation 

Let Q be a solution of the Euler or Navier-Stokes equations on the domain 
X and Q* be a discrete solution on the grid X where 


Q = 



•,j,k 


(4.4) 


and 


1 dQ 
J dt 


= R(Q)- 


Here, R is the residual and J is the transformation Jacobian 


sjiiSi o 


(4.5) 


(4.6) 


For two-dimensional thin-layer Navier-Stokes equations, the residual R can be ex- 
pressed in generalized curvilinear coordinates (£,t /) as [52], 


df d(G - G v ) 

d( dr, 

where the inviscid flux vectors F and G are 



pu 1 


pv ) 


pUu + ( X P 
pUV + (yP 

, 0.*. 

pVu + r, x P 1 
pVv + r, y P j 


( e + P)U J 


. (e + P)V j 


and U, V are the contravariant velocities defined as 


U = ZxU + t v V 


V = r, x u + T, y v. 


(4.9) 
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The viscous flux vector G v is 

f ° 1 

q 1 I 4" Vv^xy ( 

V ” J I VtTxv + Vv T vv 
[ TJxbx + Vvby . 

with shear stress, r, and heat flux, g, are written in indicial notations as 



(4.10) 


(4.11) 


and 


b - U . T + ( M °°-^ (4.12) 

b Xl + \RePr( 7 -l)J Qxi' 

where a represents the local speed of sound. The pressure is evaluated through the 
ideal gas law 

P = (7 — l)[e — p(u 2 + v 2 )/2] (4.13) 


with 7, the specific heat ratio taken to be 1.4. The molecular viscosity, /i, is calculated 
from the Sutherland’s law and Stoke’s hypothesis for bulk viscosity, A = —\p, is 
invoked. 

For a steady-state solution (i.e., t — ♦ oo)> Eq.(4.5) is reduces to 


R(Q*(P),X(P),P) = 0 


(4.14) 


where the explicit dependency of R on grid and vector of parameters P is evident. 
The parameters P control the grid X as well as the solution Q*. The fundamental 
sensitivity equation containing {^"} an( l described by Taylor et al. [13] is obtained 
by direct differentiation of Eq.(4.14) as 



(4.15) 


For a non-geometric design parameter (e.g. angle of attack or free-stream Mach 
number), the grid will not be effected by changes in P.and Eq.(4.15) reduces to 


dR 

dQ 



= 0 . 


(4.16) 
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For a geometric design parameter (e.g. camber of an airfoil), the direct dependency 
of residual on P will vanish since the changes in P would propagate through the grid. 
Equation (4.15) becomes 

ISHSMSHS)- 

It is important to notice that Eq.(4.17) is a set of linear algebraic equations , and 
the matrices [§q] and [§§;] are well understood [12]. The vector quantity {^r } 
is the solution to Eq.(4.17) given the sensitivity of the grid with respect to the 
parameters, {§jy}. A direct chain rule differentiation of {§£} results in 

{S} -[&]{*?} 

where X B designates the boundary coordinates. The vector } represents the 

boundary sensitivity which is directly related to boundary parameterization, discussed 
in the next section. It has the importance of being one of the dominant factors in cal- 
culating the sensitivity of surface forces needed for optimization process. The matrix 
l® is responsible for field grid sensitivity with respect to boundary coordinates and 
it is related to the rules which govern the grid generation algorithm. For algebraic 
generation systems, the primary components of [^^], are the interpolation functions 
which distribute the interior grid. 

An approximate version of Eq. (4.17) is suggested for predicting the steady- 
state solution changes which occur in response to small but finite boundary changes. 
The essence of using approximate analysis is to reduce the number of expensive flow 
analyses required during the optimization process. A Taylor series expansion of the 
Q*(P + AP) about P is derived as 

Q*(P + AP) w Q*(P) + AP + (4.19) 

Disregarding the higher order terms and substituting for {^} from Eq. (4.17), Eq. 


1 AQ ' 


dR 

dP 


AP 


(4.19) becomes 




(4.20) 
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where 

AQ* = Q*(P + AP) - Q*(P) 

It is essential to understand that the approximate analysis is valid as long as the 
changes in the boundary shape are small (i.e., AP 1). Figure 4.1 exhibits the non- 
linear relation between drag coefficient and thickness parameter (P = T ) for class of 
symmetrical wing-sections, enforcing the above argument. 

4.3 Surface Parameterization 

An integral part of the grid sensitivity analysis is the boundary parameter- 
ization. The process is to define the major control surfaces with independent design 
parameters. The most general parameterization of the boundaries would be to specify 
every grid point as a parameter. This conceivably could be desirable for the bound- 
aries corresponding to complex geometries to allow a design procedure to have the 
greatest possible flexibility. However, it is impractical from a computational point 
of view. Another approach would be a quasi-analytical parameterization in term of 
design variables. For instance, a class of wing-sections is specified by two camber- 
line parameters and a thickness distribution parameter; a wing is specified by several 
wing-sections; and the wing surface is interpolated from the sections. In this manner, 
an airplane component can be specified by tens of parameters instead of hundreds 
or thousands of parameters. Such physical parameterization with global boundary 
control does not possess a great deal of generality necessary for high level design and 
optimization analysis. 

A compromise between totally geometric and physical parameterization is to 
approximate the desired boundary using a spline function such as NURBS as discussed 
in Chap. 2. In the design process, using NURBS in conjunction with an interactive 
Computer Aided Design (CAD) system, would be highly advantageous [24]. After 
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prescribing an initial set of control points, the designer can pick and drag points while 
simultaneously observe the change in the shape of the surface. Although there are 
some qualifications and reservations, this approach is advocated for optimization of 
aerospace- vehicle configuration designs. 






Chapter 5 


METHOD OF SOLUTION 

5.1 Introduction 

As stated previously, the simplest way to obtain grid sensitivity is to vary 
the control parameters, one at a time, and finite difference the results. This, however, 
is proven to be computationally inefficient compared to analytical or semi-analytical 
differentiation of the grid equations. Also, the proper choice of a step size is not 
trivial and an improper choice might result in round-off error accumulation. The 
finite difference approach should only be used as the last resort when the extreme 
complexity of the grid equations dictates no other alternatives. For a less compli- 
cated grid equations, a semi-analytical approach would be more appropriate. The 
semi-analytical approach consists of analytical differentiation of the original function 
with respect to an intermediate function, the derivative of which is then evaluated 
numerically. It combines the efficiency of the analytical approach with the ease of 
implementation of the finite difference approach. 

The analytical approach to the grid sensitivity problem is evaluation of the 
grid sensitivity coefficient by direct analytical differentiation of the grid equation. For 
most cases, the grid equation is not directly differentiable, although there are schemes 
that such differentiations are feasible. The algebraic grid generation schemes, such as 
Two-Boundary Grid Generation (TBGG) used in this study, fall within that category. 
The grids are governed by explicit algebraic relations where the direct differentiation 
of the equations are not complicated. The analytical approach has the advantage of 
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being exact, thus, avoids the round-off errors associated with numerical approaches. 

There are two types of parameters involved in this study. First, there are 
the design parameters which specify the shape of the primary boundary. For classi- 
cal NACA four-digit wing-sections, the design parameters are camber, A/, location 
of camber, <7, and maximum thickness, T. For a geometric NURBS representation, 
the design parameters are the control point coordinates, D, = and their 

corresponding weights, Secondly, there are grid parameters that define other 
boundaries and the spacing between grid points. The location of far-field boundary, 
L y grid stretching parameter normal to the wall, B 3 , grid stretching parameter along 
the wall, J 9 i, and magnitude of normal derivatives along bottom and top boundaries, 
K\ and Ki, are prime examples of such parameters. The two sets of parameters are 
not functionally dependent and can be treated separately. 


5.2 Grid Sensitivity with respect to Design Parameters 

Here we express, the sensitivity of the grid with respect to the vector of de- 
sign parameters X^. As a consequence of using algebraic grid generation technique 
in which the boundary grid has the dominant effect on the interior grid, the boundary 
grid sensitivity coefficient would also be essential in influencing the interior grid sen- 
sitivity coefficient. Therefore, evaluation of the surface grid sensitivity coefficients, 
{ «&? } and 1 } , are ^e most > m P ortant P art °f analysis and directly 
dependent to the surface parameterization. 

Two distinct techniques are outlined in the Chap. 2 for wing-section param- 
eterization. The first technique, using the analytical relations of Eqs.(2.1-2.4), yields 
the classical representation of NACA four-digit wing- sections. The second technique, 
using an approximative (NURBS) relations of Eq.(2.9), would result in defining any 
free-form surface, although, the focus here would be on a wing-section geometry. For 
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practical purposes, the grid sensitivity and orthogonality at the far-field boundary 
has been ignored, (i.e., = ■ S '( s ) = °)* The evaluation of the 

grid sensitivity using analytical parameterization will be discussed first. 

A direct differentiation of Eqs.(3.9) and (3.10) with respect to Xp results in 


where 


dx 

WE 

dy 

dX D 


PS) 


g*x(r,Pj) 

dX D 


+ R(r)0l(t, Pg) 


dx\jr, P|) 
dX D 


/ 3 ?(«,P 3 ) 


gyifcZi) 

dX D 


+ R(r)0l(t, PJ) 


gyj(r,P{) 

dx D 


(5.1) 

(5.2) 


X d = {T,M,C} t . (5.3) 

The prime indicates differentiation with respect to t and can be substituted from 
Eq.(3.15). Since Xi(r, Pf) is independent of design parameters Xd, then 


a*i(r,Pl) „ 
ax D 

The x coordinate sensitivity, Eq.(5.1), can now be reduced to 


dx 

dX^ 


R{r)ft(t, PJJ) 


d(^fsin&) 

ax D ' 


Using the relation 

d , 1 du 

W" U = T7S9X D 

the x coordinate sensitivity becomes 


(5.4) 


(5.5) 


(5.6) 


dx 

dX^ 


T R{r)Pl{t,PZ)cos6 

1 + 


1 d__ 

( 9yi(^.Pf) \ 2 9Xd 

\8*i(r,pj)y 


{ 


dy ,(r,PS) l 

3xi (r, Pf ) J 


(5.7) 


The term .4- ( 9vi ( r '(*lj 1 can be evaluated by direct differentiation of Eq.(2.4). The 
SX D ( 8ii(r,Pj) J 

y coordinate sensitivity with respect to design parameters can be obtained using sim- 
ilar procedure. Equation (5.2) can be modified to 
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d f a yi (r,pjn 

T dXo \9ii(r,Pj)/ 
(5.8) 

All terms at the right hand side of Eqs.(5.7) and (5.8) can be evaluated explicitly due 
to analytical parameterization of the surface for this particular geometry. 

The grid sensitivity using the NURBS parameterization is established using 
Eqs.(5.1) and (5.2) where 




1 


dX D 


1 + 


f A 

\®*l(r.Pj ) / 


X D = (X i ,Y i ,L> i } T * = 

(5.9) 

The surface grid sensitivity coefficients, j 1 ' 1 1 are 

obtained by direct 

differentiation of Eq.(2.9) with respect to X/j as 


ax,(r, p S) dyi(r,p\) _ „ 
ax, ~ 9Y, i ’ 

(5.10) 

ai,(r, p |) aj,(r, p S) 
dVi dXi 

(5.11) 

For the weight function, we have 


dX _ 9R, lP (r)p 

doJi “ dw k 

(5.12) 

where 


x ={;:l:: p l|}- D -={?} 


and 



du k 


■,(r)N,.p(r)to>j 

E>=0 n jA t )^ 

E"=o N >. p ( r H 
k,i = 0, ...,n 


k i- * 1 
k = tj 


In the preceding equation, Ri, P (r) and N itP {r) are the basis function of the rational 
and non-rational B-splines as defined in Chap. 2. 
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5.3 Grid Sensitivity with respect to Grid Parameters 

The sensitivity of grid with respect to vector of grid parameters Xq has a 
significant impact on grid optimization and adaption. The Xq represents, as men- 
tioned earlier, the parameters that govern the surfaces other than vehicle surface and 
the parameters that control the grid spacing. Among those parameters, grid stretch- 
ing parameter along the wall B x and grid stretching normal to the wall B 3 are of 
great importance and require more attention. The grid sensitivity with respect to B\ 
is obtained as 


dx ffift Tvn 5x l( r ’ P i) | MrW(t 

w t = fit{i,Fo) + ’ l} dB x 

&y _ a0( t . R( r \B i (t p<) 9 yi( r »?i) 

TpT. ~ /W> P o)— + «\T)P dt.t',] 


(5.13) 

(5.14) 


1) , R(r W(t P^ MlL p { 
dB x ~ * u/ dB t +Wl( ’ P,) dBx 

The term can be obtained by direct differentiation of Eq.(3.11) for a bi- 

exponential distribution function as 


dxi(r, Pj) 1) e Bl (e~^^ 1)1 n < f < X 

dBx “ L (e B * — l) 2 J ' 


(5.15) 


and 

3*.(r,P?)_ f| /8B,\[ A(ef; - 1)^ - e»’(e B > A - 1) 1 x<f<l 

(5.16) 

where A = and X and Y are the inflection point coordinates for a bi-exponential 
distribution function. The quantities B 2 and can determined using New- 
ton’s iterative scheme. The term | ^ | can evaluated using the chain rule 

differentiation 

dj/i( r >Pj) [ dyi( r > p l) l [ ^x(r,Pj) ) (5 17) 

dBx ~ Uxx(r,P«)J l dBx j' 

The grid sensitivity with respect to the stretching parameter normal to wall ,B 3 , is 


dx 

m 


, w.pf.MMo) J* 

= xjfr.p,) — j- — + dl as. 


dt dB 3 
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and 


. . , phMWl+sux'i, pi' MMal d l 

+ *!(’•, Pj) g t g g 3 + «(<)*>(».?>) gt gg 3 


dy .. n.WWiL + jK^pcjSSMUiL 


(5.18) 


dB, ~ y,(r ’ P ‘) at dftj 


0i?3 


, , „ f ,a/9?((,P5) at , „ , 8«(f, PS) at 

+ »>(■•, PD—gj + $(*)»,(». P,) ^ 8ft 


(5.19) 


where 


, MMS) = 3^-4. + l 
at at 

d(3%(t, Pq) „ 2 | fif 5/?|(<,P 0 ) _ oj 2 1 
* = -61 +6i gt 31-1. 

Differentiating of Eq.(3.14) for a hyperbolic tangent stretching function yields to 

8! (t, - - l)]*o«ft(fe) - se cl. 1 (a) ia „/ l [fi.(, _ l)] ([W 

W 3 ~ 2[lanMfi‘)] 1 ( ' 

Similar developments can be extended to other grid control parameters such 
as the far-field boundary location, L , and the magnitude of orthogonality vector at 
the boundaries, K\ and K 2 . 


5.4 Flow Analysis and Boundary Conditions 

The two-dimensional thin-layer Navier-Stokes equations are solved in their 
conservation form using an upwind cell-centered finite-volume formulation. A third- 
order accurate upwind biased inviscid flux balance is used in both streamwise and 
normal directions. The finite- volume equivalent of second-order accurate central dif- 
ferences is used for viscous terms. The resulting discretization represents the residual, 
R(Q) ? at each cell depending locally on values of Q at nine neighboring cells such 
that 

R,'j(Q) — R-i,j(Qi,j > Qt'.j-l > Qi,j+1 > Qi ,j— 2) Q. j+ 2, Qi— l,j> Q«+l,i) Q«-2,j4 Qi+2,j)- 

(5.21) 
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The discretized governing equations are implicitly advanced in time using 
Euler implicit method which is unconditionally stable for all time steps according to 
Fourier stability analysis [52]. A first-order Taylor’s series expansion of the right hand 
side of Eq.(4.5) results in 


1 

JAt 


dR n 

dQ 


AQ = R"(Q) 


AQ = Q n+1 - Q n n = 1,2,3.... 


(5.22) 


where n represents the time level. The field variables for the new time level, Q 


n+l 


can be obtained by solving Eq.(5.22) for AQ. The coefficient matrix, |^||| 
square matrix with a block-banded structure of at most nine 4x4 block diagonals. An 


, is a sparse 


iterative approximate factorization (AF) algorithm have been chosen to advance the 
solution in time until 


R(Q*) « 0 


(5.23) 


where Q* are the steady-state values of the field variables. The inviscid flux vectors 
are evaluated using the upwind flux vector splitting of Van Leer [53]. For an infinitely 
large time steps, (i.e., At — > oo) the non-linear system of Eq.(5.22) may be directly 
solved using Newton’s method [14,54]. The boundary conditions are implicitly im- 
plemented within the governing equations. The airfoil surface is considered to be 
impermeable and adiabatic. A standard no-slip boundary condition with zero surface 
velocity has been selected. The pressure at the surface is evaluated using a zeroth- 
order extrapolation from the interior cells. The density is then calculated using the 
state equation, Eq.(4.13). Consequently, 


U = °, v = 0, ^=0, ^ = 0 (5.24) 

where T represents the adiabatic surface temperature and n is the unit normal vector 
of the surface. In the far-field, assuming a locally one-dimensional flow, Riemann 
invariants are employed as 


R ± =U± 


2a 

7-1' 


(5.25) 
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Adding and subtracting these two equations would yield to local normal 
velocity and speed of sound. A periodic boundary condition is applied along the wake- 
cut which resulted from a C-type grid. The effects of different boundary conditions 
on the coefficient matrix are extremely important as outline in [13,54]. For 

example, the periodic boundary condition will cause the matrix ^§qj to deviate from 
its original neat banded structure with some non-zero coefficients outside of central 
bandwidth [54]. This restructuring will greatly effect the procedure for solving the 
aerodynamic sensitivity equation to be discussed in the next section. 


5.5 Flow Sensitivity Analysis 

The flow sensitivity coefficient can now be directly obtained using 


the fundamental sensitivity equation, Eq.(4.17), as 


{£}- 


dR 


-i i 


m\ (ex 


} 


eR 

aX 


(5.26) 
, can be 


ml L3XH0P 

provided that grid sensitivity, {§p}> is known. The Jacobian matrix, 
evaluated by differentiating the discrete residuals R,j with respect to four vertices 
of each cell. Since for a typical C-type grid, the stream-wise (I-dimension) of the 
grid is significantly larger than normal (J-dimension), a J-ordering of the Eq.(5.26) 
will be implemented. This, will substantially reduce the memory requirements due 


to smaller central bandwidth of 


-l 


The quantity ^§qj can be obtained using 
a full matrix solver to account for all the non-zero contributions outside of central 
bandwidth. This, although convenient, is not practical for Navier-Stokes equations 
due to large storage requirements. An alternative would be the use of a hybrid direct 
solver with conventional relaxation strategy implemented in two steps. The coefficient 


matrix 


dR 

dQ 


can be splitted into two matrices as 

'0Rl 


dQ 


= [M] + (N) 


(5.27) 


55 


where [M] represents the banded part of the ^f|jj and [N] represents entries outside 
the band. The linear system of Eq.(5.26) can now be solved using the Richardson’s 
iteration 


{!f}‘ - + ' M >"[i] {§} - ' < 5 - 28) 

k = 1 , kmax 

where is assumed to be zero and [M] _1 can be determined by a simple con 

ventional banded matrix solver. It is evident that matrices [M] -1 , [N], and §< 
are invariant with respect to iteration index k, so they only need to be computed 
once. A single conventional under-relaxiation parameter has been used to assure the 
convergence of Eq. (5.28). For a typical design analysis of an airfoil, provides 

far more information than needed [13]. In most cases, the sensitivity of aerodynamic 
forces on the surface, such as lift and drag coefficients, are sought. The drag and lift 
coefficients are given as 

Cl = Cycosa — Cxsina _ (5.29) 

Cd = Cysina + Cxcosa (5.30) 

where a is the angle of attack. The quantities Cx a-nd Cy are the total force coeffi- 
cients along x and y directions, respectively, and can be expressed as 


Cx = ^2 C Pi {y i+ t - yi ) + C f ,(x i+ i — s<) 


i=i 


Cy = T:c p; (x,- +1 — Xi) + C/i(yi + 1 — yi) 


(5.31) 


(5.32) 


i=l 


where C p i and C/, are pressure and skin friction coefficients respectively defined as 


C P ,= 


Pi 


C], — T 


Ti 


(5.33) 


\PooUl ~ 3 ' \PooUl 

and J represents total number of boundary cells along the airfoil surface. The terms Pi 
and Tj are pressure and shear stress associated with boundary cell i and the quantity 
\pooUl\s known as dynamic pressure of the free stream. Finally, the drag and lift 


56 


sensitivity coefficients with respect to Xd 


are obtained by differentiating Eqs.(5.21) 


and (5.22) as 


8C l 

dX D 

dC D 

dX D 
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dX D 
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dX D 
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dX D 
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(5.34) 

(5.35) 


5.6 Optimization Problem 

An objective of a multidisciplinary optimization of a vehicle design is to ex- 
tremize a payoff function combining dependent parameters from several disciplines. 
Most optimization techniques require the sensitivity of the payoff function with re- 
spect to free parameters of the system. For a fixed grid and solution conditions, the 
only free parameters are the surface design parameters. Therefore, the sensitivity of 
the payoff function with respect to design parameters are needed. 

The optimization problem is based on the method of feasible directions 
[55,56] and the generalized reduced gradient method. This method has the advantage 
of progressing rapidly to a near-optimum design with only gradient information of the 
objective and constrained functions required. The problem can be defined as find- 
ing the vector of design parameters Xd, which will minimize the objective function 
/(Xd) subjected to constraints 

ffj(X D )<0 j = l,m (5.36) 

and 

X' D < X D < X U D (5.37) 

where superscripts denote the upper and lower bounds for each design parameter. 


The optimization process proceeds iteratively as 

XI = XJ- ' + 7 S" 


(5.38) 
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where n is the iteration number, S n the vector of search direction, and 7 a scalar move 
parameter. The first step is to determine a feasible search direction S n , and then 
perform a one-dimensional search in this direction to reduce the objective function 
as much as possible, subjected to the constraints. 

The present optimization strategy is based on maximizing the lift coefficient, 
Cl , in response to surface perturbation, subject to pre-determined design constraints. 
Upper and lower bounds set for each design parameter and the sensitivity derivatives 
of the objective function, and the constraint, are obtained as previously 

described. Throughout the analysis, the drag coefficient, Cd, is to be no greater than 
the value of the initial design. The strategy, illustrated in Fig. 5.1, requires that the 
grid and grid sensitivity derivatives be provided dynamically during the automated 
optimization process. 















Chapter 6 


RESULTS AND DISCUSSION 


Three test cases are considered to demonstrate the feasibility of current pro- 
cedure. For each case, the grid and flow sensitivity coefficients of the field have been 
obtained. The sensitivities of the total surface forces (i.e., Lift and Drag coefficients) 
are tabulated for optimization purposes. The first test case, a symmetrical NACA 
0012 wing-section, has been used mainly to exhibit the accuracy of the grid sensitivity 
coefficients with those obtained using the finite difference approach before proceeding 
to more challenging problems. The second test case, a NACA 8512 wing-section, has 
been used to extend the analysis to a more demanding problem involving three design 
parameters. An optimization module has been integrated into the overall procedure 
to optimize the geometry using the resultant sensitivity coefficients. Also for this 
case, some aspects of grid sensitivity with respect to grid parameters are investigated 
for grid optimization problem. The last case involves a generic representation of 
a wing-section using NURBS for surface definition. The improved design has been 
obtained employing sensitivity coefficients of the domain with respect to previously 
chosen design control points. Due to flexibility of NURBS in representing any sur- 
face, the generic wing-section test case has been used to manifest the extension of 
this approach to almost any desired geometry. 
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6.1 Case 1: NACA 0012 Wing-Section 

6.1.1 Grid Sensitivity 

The previously obtained grid, as shown in Fig.3.9, is considered for grid 
sensitivity analysis of this test case. The grid sensitivity with respect to the vector of 
design parameters Xp, are obtained using Eqs.(5.1) and (5.2) where the maximum 
thickness T is the only design parameter. 

Figure 6.1a shows the contour levels of the y-coordinate sensitivity with re- 
spect to the thickness parameter, T. The highest contour levels are, understandably, 
located in the vicinity of the chordwise location for the maximum thickness of the 
wing section. For a NACA four-digit wing section, this is positioned at about 0.3 
of the chord length from the leading edge (25). The positive and negative contour 
levels corresponding to the upper and lower surfaces are the direct consequence of 
Eq.(2.4) and the second term on the right hand side of Eq.(5.8). The sensitivity 
levels decrease when approaching the far-field boundary due to diminishing effects 
of the interpolation function /?°(f, Pq). The second term on the right hand side of 
Eq.(5.8) is responsible for the sensitivity effects due to orthogonality on the surface, 
and it is directly proportional to the magnitude of the orthogonality vector K\. The 
wake region is not sufficiently affected by any of the design parameters, and no major 
sensitivity gradient should be expected there. 

Figure 6.1b exhibits the contour levels of the x-coordinate sensitivity with 
respect to thickness parameter, T. An interesting observation can be made here 
regarding the contour levels adjacent to the surface. Unlike the y-coordinate sensi- 
tivity, the x-coordinate sensitivity is independent of design parameters at surface as 
indicated by Eq.(5.4). The contour levels resulting from Eq.(5.7) are solely due to or- 
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thogonality effects. There are some negative pockets of contour levels on the forward 
section and some corresponding positive pockets on the rear section. The dividing 
line between these pockets is located near the location of the maximum thickness (i.e., 
0.3 of the chord from the leading edge). A simple conclusion from Fig. 6.1b is that 
by increasing the thickness parameter, T, points on the forward section will move to 
the left, while at the same time, points at rear section will move to the right. 

For comparative purposes, the grid sensitivity for this case is obtained using 
the finite difference approach. The design parameters(i.e. T for this case) are per- 
turbed, one at a time, and a new grid is obtained using Eqs. (3.9) and (3.10). The 
sensitivity is then computed using a central difference approximation and the results 
are presented in Fig. 6.2 . A side by side comparison of both results indicates good 
agreement between the two approaches. 

6.1.2 Flow Sensitivity 

The second phase of the problem is obtaining the flow sensitivity coefficients 
using the previously obtained grid sensitivity coefficients. In order to achieve this, 
according to Eq.(4.17), a converged flow field solution about a fixed design point 
should be obtained. The computation is performed on a C-type grid composed of 
141 points in the streamwise direction with 101 points on the wing-section surface, 
and 31 points in the normal direction. The far-field and outer boundary were placed 
about 20 chord-length away from the airfoil. It is apparent that such a coarse grid is 
inadequate for capturing the full physics of the viscous flow over an airfoil. Therefore, 
it should be understood that the main objective here is not to produce a highly ac- 
curate flow field solution rather than to demonstrate the feasibility of the approach. 

The two-dimensional, compressible, thin-layer Navier-Stokes equations are 
solved for a free stream Mach number of M <*, = 0.8, Reynolds number Re oo = 10 6 , 
and angle of attack a = 0°. The solution is implicitly advanced in time using local 
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time stepping as a means of promoting convergence toward the steady-state. The 
residual is reduced by ten orders of magnitude as illustrated in Fig. 6.3. All compu- 
tations are performed on NASA Langley’s Cray-2 mainframe with a computation cost 
of 0.1209xl0“ 3 CPU seconds/iteration/grid point. Figures 6.4 and 6.5 demonstrate 
the pressure and Mach number contours of the converged solution with lift and drag 
coefficients of Cl = 1.53xl0~ 8 and Cd = 4.82xl0“ 2 . Due to surface curvature, the 
flow accelerates along the the upper and lower surfaces to supersonic speeds, termi- 
nated by a weak shock wave behind which it becomes subsonic. 

The sensitivity coefficient, is obtained by previously described iter- 

ative strategy of Eq.(5.28). The average error has been reduced by three orders of 
magnitude and the convergence history is shown in Fig. 6.6. The sensitivities of the 
aerodynamic forces, such as drag and lift coefficients with respect to thickness pa- 
rameter T, are obtained utilizing Eqs. (5.29-5.35) with results presented in Table 6.1. 
Again, for comparison , a finite difference approximation has been implemented to 
validate the results. A nominal perturbation of 10 -4 for design parameter T has been 
chosen and the corresponding results are included in Table 6.1. The good agreement 
between the two sets of numbers verifies the accuracy of the approach. It is apparent 
from Table 6.1 that while drag is extremely sensitive to the changes in wing-section 
thickness, the lift is almost insensitive, due to symmetrical nature of flow for this 
case. 


Table 6.1 Lift and Drag sensitivities with respect to design parameter T 


NACA 0012 

Direct Approach 

Finite Difference 

Design Parameters 

oUl 

dXn 

dCp 

dXn 

dCt. 

dXn _ 

dUn 

dXn 

Maximum Thickness 

—2.68x10“® 

0.709 

-4.75x10“® 

0.707 






(a) Y-Coordinatc 


0435 * 1 * 
0.373557 
0311207 
0.249030 
0.100770 
0.124510 
0.062259 
7.450E 9 
-0.062259 
-0.124510 
-0 186770 
-0240030 
0.311297 
-0.373557 
-0 435010 



0025997 
0017061 
0009924 
0.001087 
-0006140 
-0.014105 
-0022222 
-0 030258 
-0038295 
-0.046332 
-0.054360 
-0062405 
-0 070441 
-0070470 
-0 006515 


(b) X-Coordinate 

Fig. 6.2 Coordinate sensitivity with respect to maximum thickness T (FD) 








Fig. 6.5 Mach number contours for NACA 0012 wing-section 
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6.2 Case 2: NACA 8512 Wing-Section 

0.2.1 Grid Sensitivity 

The second test case considered is the NACA 8512 cambered airfoil. Again, 
the previously obtained grid, as shown in Fig. 3.10, is considered for grid sensitivity 
analysis. Figure 6.7 shows the coordinate sensitivity with respect to parameter T with 
characteristics similar to the previous symmetrical airfoil case. Figure 6.8a represents 
the y-coordinate sensitivity with respect to camber, Af. It appears that the highest 
sensitivity contour levels are located at the chordwise location of camber, C (i.e., 0.5 
of the chord-length). The contour levels decrease toward the far- field boundary, again 
as a consequence of interpolation function. However, unlike Fig. 6.7a, they possess 
positive values on both upper and lower surfaces. Consequently, an increase in cam- 
ber, M, shifts all the points upward. Again, minimum activity can be detected in the 
wake region. Figure 6.8b shows the x-coordinate sensitivity contours with respect to 
camber, M. Here, as in Fig. 6.7b, the sensitivities are minimum on the surface of 
the wing-section. There is a small gradient on the forward section, but by far the 
strongest gradient is in the rearward section due to orthogonality effects. 

Figure 6.9a illustrates the y-coordinate sensitivity with respect to camber 
location, C. A dividing line between positive and negative contour levels appears 
near the chordwise position of the camber. Like previous cases, there is no significant 
activity in the wake region. The result indicates that a positive change of C will cause 
the movement of points downward on the forward section, while at the same time, the 
points on the rear section will respond by moving upward. Figure 6.9b illustrates the 
x-coordinate sensitivity with respect to camber location C. The two major features 
are attributed to chordwise location of the camber and the orthogonality effects on 
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the tail section. It is interesting to notice that the sensitivity level for camber location 
is considerably less than the other two design parameters. 

Another important aspect of grid sensitivity analysis is to investigate the 
sensitivity of a grid with respect to grid parameters. Such information would be most 
helpful in grid optimization, therefore, improving the flow solution and overall opti- 
mization process. Figure 6.10 illustrates the grid sensitivity with respect to stretching 
parameter, £? 3 , where the hyperbolic behavior of contour levels are apparent. The 
sensitivity of grid with respect to grid distribution parameter around an airfoil, Bi, is 
shown in Fig. 6.11. The y-coordinate sensitivity does not register any large gradient 
field except on the front. The x-coordinate, being the distribution axes, do however 
show some interesting contour levels on the surface. The sensitivities with respect 
to orthogonality parameter at the surface, are presented on Fig. 6.12. Figure 
6.13 exhibits the grid sensitivity of the domain with respect to the outer boundary 
location, L. The effect of outer boundary location on the solution is a subject of 
extensive research and deserves a more comprehensive investigation [33]. 

0.2.2 Flow Sensitivity and Optimization 

Using free stream conditions of = 0.7 , Re = 10 6 , and a = 0°, a 
converged flow field solution is obtained. As in the previous case, a C-type grid of 
141x31 is used and the far-field boundary is placed about twenty chord-length away 
from the surface. The residual is reduced by ten orders of magnitude using 3000 
implicit Euler time iterations. Figures 6.14 and 6.15 illustrate the resulting pressure 
and Mach number contours. As a consequence of camber, flow accelerates along the 
upper surface to supersonic speeds. As the flow travels further, it encounters an ad- 
verse pressure gradient due to decreasing curvature, resulting in a weak shock and 
subsonic flow. Figure 6.16 shows the surface pressure coefficient C p , where lift and 
drag coefficients are Cl = 0.611, and Cd = 0.094. 
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The surface aerodynamic sensitivity coefficients with respect to the vector of 
design parameters Xr>, are obtained and presented in Table 6.2. Figure 6.17 demon- 
strates the convergence behavior of Eq. (5.28) for P = {T,M,C} T . It is apparent 
that a combination of increasing camber and decreasing thickness would greatly im- 
prove the Lift/Drag ratio for this particular problem. The effects of camber location 
C is almost negligible and might be ignored in optimization process in order to de- 
crease the computational costs. A comparison with finite difference approximation s 
reveals the accuracy of the this direct approach. It is imperative to understand that 
due to the non-linear nature of Navier-Stokes equations, the results of Table 6.2 is 
only valid for current specified wing-section and flow condition. Extrapolating from 
these coefficients to predict flow solutions for different families of wing-sections would 
be greatly erroneous and unreliable. Such predictions are only valid as long as the 
design perturbations are in the same order as the finite-difference step size used in 
Table 6.2. 

Table 6.3 displays similar results for the vector of grid parameters, Xq. 
The far-field boundary location L has the greatest effect on lift and drag, followed 
by surface grid orthogonality parameter, K\. An underlying effort in generating a 
suitable grid for flow analysis is to minimize the gridding effects on the solution. 
Inspecting Table 6.3 reveals that the solution is apparently grid sensitive; therefore, 
not particularly suitable for this geometry. This grid dependency of solution may be 
mainly blamed on the coarseness of the grid and the location of far-field boundary. 
Infact, according to Table 6.3, the far-field boundary should be placed further away 
from wing-section in order to achieve better Lift/Drag ratio. This, will undoubtedly 
contribute to the coarseness of the grid, and ultimately to solution instability. In- 
creasing the grid density, although may alleviate the problem, but will cause a surge 
in computational costs. Another improvement may be to modify the orthogonality 
vector in response to high slope regions of the surface. Also, a modest increase in 
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stretching parameter might produce better results as indicated in Table 6.3. The 
grid distribution on the surface appears to be sufficient due to small value of its sen- 
sitivity coefficients. An ideal situation would be to alter the grid parameters until 
the sensitivity coefficients of Table 6.3 all approach to zero. A combination of the 
above recommendations although might produce such grid, but would be tedious and 
time consuming. An alternative is to optimize each of these grid parameters using 
similar optimization techniques as devised for design parameters. In this respect, an 
optimum grid can be constructed for this particular example which also can be used 
in design optimization cycles. Since the main objective here is to demonstrate the 
validity of grid sensitivity module and its integration into the optimization loop, these 
deficiencies will be overlooked for the present time. 

The design optimization strategy of Fig. 5.1 is applied to three design pa- 
rameters of T, M , and C. The intention is to maximize the Lift/Drag ratio subjected 
to appropriate boundary perturbation. The upper and lower bounds for design pa- 
rameters are assigned as 

0.08 < T < 0.16, 0.04 <M< 0.12, 0.3 < C < 0.7. (6.1) 

The computation is performed dynamically using the design optimization 
module of ADS (Automated Design Synthesis) as outlined in [56]. The optimum 
design is achieved after 11 non-linear optimization cycles and a total of 6800 CPU 
seconds on Cray-2 mainframe. A simple break down of computational costs indicates 
that about 80% of this processing time is spent on flow analyses , as is the case for most 
flow solvers. Table 6.4 shows the improvement in Lift/Drag ratio for this particular 
wing-section. The corresponding design parameters and their optimum values are 
included in Table 6.5. As expected from sensitivity coefficients, maximum thickness 
T, and maximum camber M , had the greatest effect on this optimization process. 
Figure 6.18 illustrates the comparison between initial and optimized wing-section. 
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Table 6.2 Lift and Drag sensitivities with respect to vector design parameter X D 


NACA 8512 

Direct A 

jproach 

Finite Difference 

mi- r-" 

Design Parameters 

dV l 
dXn 

dC'n 

Mu 

dXn 

oo n 
dXn 

Maximum Thickness T 

-4.23 

0.442 

-4.23 

0.442 

Maximum Camber M 

+7.84 

0.322 

+7.83 

0.322 

Maximum Camber location C 

-8.34xl0 -5 

-4.3xl0- 3 

-8.5xl0 -3 

-4.1xl0" 3 


Table 6.3 Lift and Drag sensitivities with respect to vector of grid parameters X G 


NACA 8512 

Direct Approach 

Finite Difference 

Grid Parameters 

dXn 

dCn 

Ma - 

dXa 

d -H- 

Grid Stretching 

l.lxlO" 2 

1.5xl0 -3 

1.08xl0" 2 

1.5xl0" 5 

Grid Distribution 

— 2.2xl0 -3 

— 8.6xl0 -4 

— 2.2xl0 -3 

-8.58xl0" 4 

Surface Orthogonality 

+3.2xl0 -2 

8.13xl0 -3 



Far-held Boundary Location 

USSSSM 





Table 6.4 Comparison of initial and optimized performance variables 


Performance 

mmm 

Optimum 

Percent 

Variables 

BM^I 

Design 

Change 

C L 1 

0.611 

0.746 

+22 

Cd 

niiniHi 

0.089 

-5.3 

Lift/Drag Ratio 


8.38 

+28.9 


Table 6.5 Comparison of initial and optimized design parameters 


Design 

Initial 

Optimum 

Percent 

Parameters 

Design 

Design 

Change 

Max. Thickness T 

0.12 


-33.33 

Max. Camber M 

0.08 


+22.5 

Location of Camber C 

0.5 

0.55 

+10 
































73 



0 034711 
0 060729 
0 004677 
0 730624 
0 674572 
0 60052 
0544467 
0 470415 
0 414363 
0340311 
0 264256 
0 210206 
0 154154 
0 060101 
0 024040 



F 

E 

D 

C 

B 

A 

1 


0 005107 
0061223 
0 067330 
0 053455 
0 030570 
0 025666 
0011002 
0002061 
0 015966 
0 029650 
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Fig. 6.8 Coordinate sensitivity with respect to maximum camber M. 
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Fig. 6.11 Coordinate sensitivity with respect to grid distribution parameter B i 
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Fig. 6.12 Coordinate sensitivity with respect to grid orthogonality, parameter K\ 
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6.3 Case 3: Generic Wing-Section 

0.3.1 Grid Sensitivity 

The grid sensitivity with respect to design parameters using the NURBS 
parameterization is discussed in this section. A generic wing-section is devised by 
employing Eq.(2.9) and seven pre-specified control points and weights (option 3). 
Figure 6.19 represents the wing-section and the control polygon using cubic basis 
functions of Fig. 2.9. The control points are numbered counter-clockwise, starting 
and ending with control points (0 and 6), assigned to the tail of the wing-section. As a 
consequence of Eq.(5.9), the total number of design parameters now jumps to 21 (i.e., 
three design parameters/control point). Depending on desired accuracy and degree of 
freedom for optimization, the number of design parameters could be reduced for each 
particular problem. For the present case, such reduction is achieved by considering 
fixed weights and chord-length. Out of the remaining four control points with two 
degrees of freedom for each, control points 1 and 5 would probably have the greatest 
impact due to their camber-like positions. The number of design parameters now 
reduced to four with X D - {X u Yi,X 6 ,Y s } t with initial values specified in Fig. 6.19. 
In accordance with Eq.(5.10), the non-zero contribution to the surface grid sensitivity 
coefficients of these control points, are the basis functions f?i, 3 (r) and f? 5 , 3 (r), shown 
in Fig. 6.20. Figure 6.21 illustrates the field-grid sensitivity with respect to design 
parameter Y\. It is interesting to notice the similarities between contour patterns 
of Fig. 6.21a and 6.8a with one obvious difference. The sensitivity gradients are 
restricted only to the region influenced by the elected control point. This locality fea- 
ture of the NURBS parameterization makes it a desirable tool for complex design and 
optimization when sometimes only a local perturbation of the geometry is warranted. 


0 - 2 - 
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Figure 6.22 represents similar results for design control point 5 where the sensitivity 
gradients are restricted to the lower portion of domain. 

0.3.2 Flow Sensitivity and Optimization 

The previous free-stream conditions of case 2 has been used to advance the 
solution in time. The far-held boundary is again placed at 20 chord-length away with 
specified boundary conditions outlined in section 5.4. The previously defined C-type 
grid of 141x31 grid points is used and the residual is reduced by ten orders of magni- 
tude. Flow characteristics similar to case 2 are detected with lift and drag coefficients 
of Ci = 0.402 and Co = 0.063. The lift and drag sensitivities with respect to Xp 
are presented in Table 6.6. The finite -difference comparison of sensitivity coefficients 
has been avoided due to already verified sensitivity module. An inspection of Table 
6.6 indicates the substantial influence of parameters Y\ and Y$ on the aerodynamic 
forces acting on the surface. The upper and lower bounds for these design parameters 
are assigned as 

0.2 < Xi < 0.7, -0.1 < Y\ < 0.5, 0.2 < X s < 0.7, -0.1 < Y s < 0.2. (6.2) 

The optimum design is achieved after 17 optimization cycles and total of 
8807 Cray-2 CPU seconds. Comparing with case 2, an almost 30% CPU increase can 
be attributed to the addition of extra design parameter. These computational cost 
overruns make minimizing the number of design parameters in CFD optimization 
essential. Table 6.7 highlights the initial and final values of lift and drag coefficients 
with a 208% improvement in their ratio. Table 6.8 represents the initial and optimum 
design parameters with parameters Y\ and Y& having the largest change as expected. 
The history of design parameters deformation during the optimization cycles appears 
in Fig. 6.23, where the oscillatory nature of design perturbations during the early 
cycles are clearly visible. Figure 6.24 compares the initial and optimum geometry of 
the wing-section. 
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Several observations should be made at this point. First, although control 
points 1 and 5 demonstrated to have substantial influence on the design of the wing- 
section, they are not the only control points affecting the design. Infact, control points 
2 and 4 near the nose might have greater affect due to sensitive nature of lift and drag 
forces on this region. The choice of control points 1 and 5 here was largely based on 
their camber like behavior. A complete design and optimization should include all the 
relevant control points (e.g., control points 1, 2, 4, and 5). For geometries with large 
number of control points, in order to contain the computational costs within a reason- 
able range, a criteria for selecting the most influential control points for optimization 
purposes should be implemented. This decision could be based on the already known 
sensitivity coefficients, where control points having the largest coefficients could be 
chosen as design parameters. Secondly, the optimum wing-section of Fig. 6.24 is only 
valid for this particular example and design range. As a direct consequence of the 
non-linear nature of governing equations and their sensitivity coefficients, the validity 
of this optimum design would be restricted to a very small range of the original design 
parameters. The best estimate for this range would be the finite-difference step size 
used to validate the sensitivity coefficients (i.e., 10” 3 or less). All the wing-sections 
with the original control points within this range should conform to the optimum 
design of the Fig. 6.24, while keeping the grid and flow conditions fixed. 
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Table 6.6 Lift and Drag sensitivities with respect to vector of grid parameters Xp 


Generic Wing-Section 

Direct Approach 

Design Parameters: 

oXn 

w- 

°xq 


-0.297 

IffnfnT'lM 

Vi — I 

-5.107 

0.549 I 

x 5 

0.15 


r 8 

2.609 

0.287 | 


Table 6.7 Comparison of initial and optimized performance variables 


Performance 

Variables 

Initial 

Design 

Optimum 

Design 

Percent 

Change 

C L 

0.402 

0.845 

+110.1 

Cd 

0.063 

0.043 

-31.7 

Lift/Drag Ratio 

6.38 

19.65 

+208 


Table 6.8 Comparison of initial and optimized design parameters 


Design 

Parameters 

Xx 

Yx 

*6 

Y s 


Initial 

Design 

05 

0.2 

0.5 

0.05 


Optimum 

Design 

0374 

0.134 

0.414 

0.069 


Percent 

Change 

-25.2 

-33 

-17.2 

+38 
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Fig. 6 .22 Coordinate sensitivity with respect to Y s . 








Chapter 7 


CONCLUSIONS AND RECOMMENDATIONS 


An algorithm is developed to obtain the grid sensitivity with respect to de- 
sign parameters for aerodynamic optimization. The algebraic Two-Boundary Grid 
Generation (TBGG) scheme has been directly differentiated with respect to design 
parameters. This formulation has the benefits of being exact, efficient, and inexpen- 
sive. The methodology is applied first to a symmetrical wing-section where the grid 
sensitivity coefficients have been validated by comparing well with finite difference 
approach (case 1). The wing-section geometry is defined either analytically using a 
combination of camber and thickness distributions, or geometrically using the NURBS 
approximation of the surface. The next test case involving a cambered wing-section 
is devised to couple the present grid sensitivity module with the newly developed 
flow analysis and sensitivity module of Taylor et al. [13], case 2. The aerodynamic 
sensitivity coefficients again compared well with finite difference results verifying the 
accuracy of grid sensitivity coefficients and their flow counterparts. Another applica- 
tion of this scheme, grid sensitivity with respect to grid parameters, has been obtained 
for grid-optimization purposes. The algorithm is then used for the main application, 
which is to optimize a generic wing-section using geometric NURBS parameterization 
of the surface (case 3). A substantial increase in aerodynamic performance variables 
enforces the feasibility of this approach for high level design and optimization. 

It is evident that grid sensitivity plays a significant role in the aerodynamic 
optimization process. The algebraic grid generation scheme presented here intended 
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to demonstrate the elements involved in obtaining the grid sensitivity from an alge- 
braic grid generation system. Each grid generation formulation requires considerable 
analytical differentiation with respect to parameters which control the boundaries as 
well as the interior grid. It is implied that aerodynamic surfaces, such as a wing- 
section considered here, should be parameterized in terms of design parameters. Due 
to the high cost of aerodynamic optimization process, it is imperative to keep the num- 
ber of design parameters as low as possible. Analytical parameterization, although 
facilitates this notion, has the disadvantage of being restricted to simple geometries. 
A geometric parameterization such as NURBS, with local sensitivity, has been advo- 
cated for more complex geometries. 

Future investigations should include the implementation of present approach 
using larger grid dimensions, adequate to resolve full physics of viscous flow analysis. 
A grid optimization mechanism based on grid sensitivity coefficients with respect to 
grid parameters should be included in the overall optimization process. An optimized 
grid applied to present geometry, should increase the quality and convergence rate 
of flow analysis within optimization cycles. Other directions could be establishing a 
link between geometric design parameters (e.g., control points and weights) and basic 
physical design parameters (e.g., camber, thickness). This would provide a consistent 
model throughout the analysis which could easily be modified for optimization. Also, 
the effects of including all the relevant control points on the design cycles should be 
investigated. Another contribution would be the extension of the current algorithm 
to three-dimensional space for complex applications. For three-dimensional appli- 
cations, even a geometric parameterization of a complete aerodynamic surface can 
require a large number of parameters for its definition. A hybrid approach can be 
selected when certain sections or skeletal parts of a surface are specified analytically 
or with NURBS and interpolation formulas are used for intermediate surfaces. 
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APPENDIX A 


TRANSFINITE INTERPOLATION WITH 
L AGRAN GIAN BLENDING FUNCTION 


A.l Surface Grid Generation 


The dual-block grid topology of Fig. 3.3 has been selected and applied to a 
generic airplane configuration. The next step is to generate the solid surfaces either 
by a set of analytical functions or by a set of cross-sections. The former requires 
no interpolation while the latter requires some sort of bi-directional interpolation. 
A powerful interpolation choice would be the use of the NURBS function for inter- 
mediate surfaces. For this geometry, a fuselage with circular cross-section has been 
devised using polar coordinates. The wing and tails are derived by using NApA 


four-digit wing-sections as a selection of cross-sections and then interpolate for inter- 
mediate surfaces. Figure A.l shows the corresponding surface grid. For remaining 
non-physical surfaces, the best approach is to decompose the region to a number of 
sub-regions as illustrated in Fig. A. 2. Although this sub-division is arbitrary, it is 
a good idea to sub-divide along computational coordinates. The grid for each sub- 
region, F{£,t)) = {x{(,r]),y{(,ri)} T , can be computed using a two-step Transfinite 
Interpolation as 

L p 
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where the /4” and J3" are the known coordinate lines on the surface with their deriva- 
tives, 
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and <4 n) (£) and /3^\t]) are the univariant blending functions, 
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These functions are subjected to the following conditions 
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These conditions permit to reproduce the input boundaries. 

Selection of the blending function depends on the number of specified bound- 
aries. One choice for the blending function is the Lagrangian interpolation which 
satisfies the Cardinal conditions. For example, if some lines in the £- direction are 
given at £1, £2* - • • * then the blending function a can be defined as, 


(£\ _ TT (£ £>) 

3*1 
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If only two boundaries are defined in one computational direction, then the La- 
grangian interpolation would convert to a simple linear interpolation 


«.«) = 


(6 - 0 
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This works if the boundaries do not contain sharp discontinuities. Otherwise, these 
discontinuities will propagate into the interior regions. One way to solve this problem 
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wa. 


is to construct a blending function that has very small value away from the boundaries. 
For example, the following blending functions have these criteria. 

- 1 


Q i(0 = 




MO = 


{A. 11) 


e K _ l ’ — t K - 1 

where K is a negative number greater than one. The larger the K is the lesser the 
discontinuity will propagate. A similar blending function can be constructed for the r/ 
direction. Some results of this formulation applied to symmetry and outer boundary 
planes are shown in Figs. A. 3 and A. 4. 


A. 2 Volume Grid Generation 


In general, decomposition of the physical domain produces several blocks. 
Each block is usually defined by six sides, and each side can be defined by either a 
surface, plane, line, or a point. If one side of a block collapses to a line or a point, 
then there would be a singularity in the block. In some instances, a block may have 
been defined by more or less than six surfaces. Once the surface are defined, the 
interior grid can be computed by any standard grid generation technique. In this 
study, an oscillatory- transfinite interpolation has been used to generate the interior 
grid points. 

Once the boundary surfaces are known, then it is possible to 

generate the interior grid by transfinite interpolations. In a general form, the transfi- 
nite interpolation (or uni variant interpolation) can be expressed by a vector F(OViO 
as 


= { x {(tV^0^ y(0 •?>()> *(&>?.<)} ’ 
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where the A?, B? and C? are the known surface location and their derivatives. Figure 
A. 5 displays the interior grid for a constant-I surface. 











